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[57] ABSTRACT 

A method for imaging objects in highly scattering turbid 
media. According to one embodiment of the invention, the 
method involves using a plurality of intersecting source/ 
detectors sets and time-resolving equipment to generate a 
plurality of time-resolved intensity curves for the diffusive 
component of light emergent from the medium. For each of 
the curves, the intensities at a plurality of times are then 
inputted into the following inverse reconstruction algorithm 
to form an image of the medium: 

X {k+l)T = [y t W + X {k)T A ][W t W + A ] -1 

wherein W is a matrix relating output at source and detector 
positions r s and r d , at time t, to position r, A is a regular- 
ization matrix, chosen for convenience to be diagonal, but 
selected in a way related to the ratio of the noise, <nn>to 
fluctuations in the absorption (or diffusion) X- that we are 
trying to determine: 

A—Ajbij with Ay=<nn>/<AXjAXj> 

Y is the data collected at the detectors, and X* is the kth 
iterate toward the desired absorption information. An 
algorithm, which combines a two dimensional (2D) matrix 
inversion with a one-dimensional (ID) Fourier transform 
inversion is used to obtain images of three dimensional 
hidden objects in turbid scattering media. 

51 Claims, 13 Drawing Sheets 
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TIME-RESOLVED DIFFUSION 
TOMOGRAPHIC 2D AND 3D IMAGING IN 
HIGHLY SCATTERING TURBID MEDIA 

CROSS-REFERENCE TO RELATED 
APPLICATIONS 

The present application is a continuation-in-part of a 
pending U.S. patent application Ser. No. 08/797,028, filed 
Feb. 7, 1997 in the names of Robert R. Alfano, Wei Cai, 
Feng Liu, Melvin Lax and Bidyut B. Das and entitled 
TIME-RESOLVED DIFFUSION TOMOGRAPHIC 2D 
AND 3D IMAGING IN HIGHLY SCATTERING TURBID 
MEDIA, which in turn is a continuation-in-part of pending 
U.S. patent application Ser. No. 08/618,471, filed Mar. 18, 
1996, both of which is incorporated herein by reference. 

STATEMENT AS TO RIGHTS TO INVENTIONS 
MADE UNDER FEDERALLY-SPONSORED 
RESEARCH AND DEVELOPMENT 

This invention was made, in part, with Government 
support awarded by NASA and ONR. The Government may 
have certain rights in this invention. 

BACKGROUND OF THE INVENTION 

The present invention relates generally to the imaging of 
objects in highly scattering turbid media and more particu- 
larly to a novel technique for imaging objects in highly 
scattering turbid media. 

As can readily be appreciated, there are many situations 
in which the detection of an object present in a highly 
scattering turbid medium is highly desirable. For instance, 
the detection of a tumor embedded within a tissue is one 
such example. Although X-ray techniques do provide some 
measure of success in detecting objects in turbid media, they 
are not well-suited for detecting very small objects, e.g., 
tumors less than 1 mm in size, or for detecting objects in 
thick media. In addition, X-ray radiation can present safety 
hazards to a person exposed thereto. 

An alternative technique used to detect objects in turbid 
media is transillumination. In transillumination, visible or 
near infrared (NIR) light is incident on one side of a medium 
and the light emergent from the opposite side of the medium 
is used to form an image. Objects embedded in the medium 
typically absorb the incident light and appear in the image as 
shadows. Unfortunately, the usefulness of transillumination 
as a detection technique is severely limited in those 
instances in which the medium is thick or the object is very 
small. This is because light scattering within the medium 
contributes to noise and reduces the intensity of the unscat- 
tered light used to form the image shadow. 

To improve the detectability of small objects located in a 
turbid medium using transillumination, many investigators 
have attempted to selectively use only certain components of 
the transilluminating light signal. This may be done by 
exploiting the properties of photon migration through a 
scattering medium. Photons migrating through a turbid 
medium have traditionally been categorized into three major 
signal components: (1) the ballistic (coherent) photons 
which arrive first by traveling over the shortest, most direct 
path; (2) the snake (quasi-coherent) photons which arrive 
within the first St after the ballistic photons and which 
deviate, only to a very slight extent, off a straight-line 
propagation path; and (3) the diffusive (incoherent) photons 
which experience comparatively more scattering than do 
ballistic and snake photons and, therefore, deviate more 
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considerably from the straight-line propagation path fol- 
lowed by ballistic and snake photons. 

Because it has been believed that ballistic and snake 
photons contain the least distorted image information and 
5 that diffusive photons lose most of the image information, 
efforts to make transillumination work most effectively with 
turbid media have focused on techniques which permit the 
selective detection of ballistic and snake photons while 
rejecting diffusive photons. This process of selection and 
io rejection has been implemented in various time-gating, 
space-gating and time/space -gating techniques. Patents, 
patent applications and publications which disclose certain 
of these techniques include U.S. Pat. No. 5,140,463, inven- 
tors Yoo et al., which issued Aug. 18, 1992; U.S. Pat. No. 
15 5,143,372, inventors Alfano et al., which issued Aug. 25, 
1992; U.S. Pat. No. 5,227,912, inventors Ho et al., which 
issued Jul. 13, 1993; U.S. Pat. No. 5,371,368, inventors 
Alfano et al., issued Dec. 6, 1994; Alfano et al., “Photons for 
prompt tumor detection,” Physics World, pp. 37-40 (January 
20 1992); Wang et al., “Ballistic 2-D Imaging Through Scat- 
tering Walls Using an Ultrafast Optical Kerr Gate,” Science, 
Vol. 253, pp. 769-771 (Aug. 16, 1991); Wang et al., “Kerr- 
Fourier imaging of hidden objects in thick turbid media,” 
Optics Letters, Vol. 18, No. 3, pp. 241-243 (Feb. 1, 1993); 
25 Yoo et al., “Time-resolved coherent and incoherent compo- 
nents of forward light scattering in random media,” Optics 
Letters, Vol. 15, No. 6 , pp. 320-322 (Mar. 15, 1990); Das et 
al., “Ultrafast time-gated imaging in thick tissues: a step 
toward optical mammography,” Optics Letters, 18(13) 
30 : 1092-4 (1993); Chen et al., “Two-dimensional imaging 
through diffusing media using 150-fs gated electronic holog- 
raphy techniques,” Optics Letters, Vol. 16, No. 7, pp. 
487-489 (Apr. 1, 1991); Duncan et al., “Time-gated imaging 
through scattering media using stimulated Raman 
35 amplification,” Optics Letters, Vol. 16, No. 23, pp. 
1868-1870 (Dec. 1, 1991), all of which are incorporated 
herein by reference. 

Of the above-listed art, Wang et al., “Kerr-Fourier imag- 
ing of hidden objects in thick turbid media,” Optics Letters, 
40 Vol. 18, No. 3, pp. 241-243 (Feb. 1, 1993) is illustrative of 
transillumination techniques which selectively use the bal- 
listic and/or snake components of light. In this article, there 
is disclosed a time/space-gating system for use in imaging 
opaque test bars hidden inside a 5.5 cm-thick 2.5% Intralipid 
45 solution. The disclosed system includes three main parts: a 
laser source, an optical Kerr gate and a detector. The laser 
source is a picosecond mode-locked laser system, which 
emits a 1054 nm, 8 ps laser pulse train as the illumination 
source. The second harmonic of the pulse train, which is 
50 generated by transmission through a potassium dihydrate 
phosphate (KDP) crystal, is used as the gating source. The 
illumination source is sent through a variable time-delay and 
is then used to transilluminate, from one side, the turbid 
medium containing the opaque object. The signal from the 
55 turbid medium located at the front focal plane of a lens is 
collected and transformed to a Kerr cell located at its back 
focal plane (i.e., the Fourier-transform spectral plane of a 4F 
system). That portion of the Kerr cell located at the focal 
point of the 4F system is gated at the appropriate time using 
60 the gating source so that only the ballistic and snake com- 
ponents are permitted to pass therethrough. The spatial- 
filtered and temporal-segmented signal is then imaged by a 
second lens onto a CCD camera. 

Although techniques of the type described above, which 
65 selectively use ballistic and snake photons to image objects 
in turbid media, have enjoyed a modicum of success, such 
techniques have been limited by the fact that detected light 
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signals derived from ballistic and snake photons are typi- 
cally rather weak, due to the proportionately small number 
of transilluminated ballistic and snake photons. This prob- 
lem is further exacerbated in those instances in which the 
turbid medium is thick and the likelihood of substantial 
scattering increases. 

Accordingly, because diffusive photons constitute the 
greatest component of the transilluminated light signal, it 
would be highly desirable to make use of the diffusive 
component of the light signal in forming an image via 
transillumination. This objective is made difficult, however, 
by the fact that diffusive photons tend to traverse a medium 
along ill-defined paths. One approach to this problem has 
been to invert the experimental scattering data obtained from 
various points in the medium using some inverse algorithm 
and reconstruction approach. This approach is often called 
diffusion tomography since diffusion or scattering is the 
dominant factor in the problem. In diffusion tomography, 
one produces an internal map of the scattering medium using 
the scattered signals and a mathematical inversion algo- 
rithm. The inversion is based upon the physical and math- 
ematical principles governing photon propagation in turbid 
media. Both time-resolved data and frequency domain data 
can be used for reconstruction. Examples of diffusion 
tomography techniques include Arridge, “The Forward and 
Inverse Problems in Time Resolved Infra-Red Imaging,” 
Medical Optical Tomography: Functional Imaging and 
Monitoring SPIE Institutes, Vol. IS11, G. Muller ed., 31-64 
(1993); Singer et al., “Image Reconstruction of the Interior 
of Bodies That Diffuse Radiation,” Science , 248:990-3 
(1993); Barbour et al., “A Perturbation Approach for Optical 
Diffusion Tomography Using Continuous-Wave and Time- 
Resolved Data,” Medical Optical Tomography: Functional 
Imaging and Monitoring SPIE Institutes, Vol. IS11, G. 
Muller ed., 87-120 (1993); M. Patterson et al., SPIE, 1767, 
372 (1992); J. Schotland et al., App. Opt., 32, 448 (1993), all 
of which are incorporated herein by reference. 

The foregoing diffusion tomography techniques do not 
lead to a resolution that is better than about 5-10 mm. 
Moreover, these techniques are time-consuming and do not 
readily lend themselves to real-time use. 

SUMMARY OF THE INVENTION 

It is an object of the present invention to provide a novel 
technique for imaging objects present in highly scattering 
turbid media. 

It is another object of the present invention to provide a 
technique as described above that overcomes at least some 
of the disadvantages associated with existing techniques for 
imaging objects in highly scattering turbid media. 

Additional objects, as well as features and advantages, of 
the present invention will be set forth in part in the descrip- 
tion which follows, and in part will be obvious from the 
description or may be learned by practice of the invention. 

The present invention is based on the present inventors’ 
discovery that objects hidden in highly scattering turbid 
media can be imaged by a novel, fast and noise-resistant 
inverse method which involves illuminating a highly scat- 
tering turbid medium with a pulse of light, the light emergent 
from the highly scattering turbid medium consisting of a 
ballistic component, a snake-like component and a diffusive 
component; determining the intensity of said diffusive com- 
ponent at a plurality of points in time; and using said 
intensity determinations to form an image of the object in the 
highly scattering turbid medium. 

The aforementioned illuminating and determining steps 
preferably include the use of a light source and one or more 
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detectors, said detectors preferably including means, such as 
a streak camera, for temporally resolving the diffusive 
component of the emergent light from the highly scattering 
turbid medium. More preferably, two or more source/ 
5 detector sets are arranged in an intersecting manner to 
enable a 2-D or 3-D image to be formed. 

The aforementioned illuminating and determining steps 
more preferably include the use of one or more two- 
dimensional area detectors, such as a charge-coupled device 
10 (CCD) or an infrared area camera. The detectors preferably 
include means to amplify the emergent light pulse from the 
turbid medium and to carve out any desired temporal slice 
of the pulse. The output of the detector is a shadow- image 
or a two-dimensional time -dependent intensity distribution, 
15 I(x,y,t), where the direction of incidence of light pulse is 
taken to be along the z direction. One such detector con- 
figuration comprises a gated image intensifier coupled to a 
CCD camera. This particular detection and imaging scheme 
is a specific embodiment of the electronic time-gating tech- 
20 nique and will hereinafter be referred to as an electronic 
time-gated imaging camera (TGIC) system. The duration 
and the position of the gate is variable from a few picosec- 
onds to a few nanoseconds. Such a system presently used by 
the inventors provides a minimum gate duration of 80 ps. 
25 The gate position can be varied in 25 ps steps. 

Preferably, the following inverse reconstruction algorithm 
is used to form an image of the turbid medium using the 
aforementioned intensity determinations: 

30 T T - 1 

X {k+ 1) =[Y T W + X {k) A][W t W + A] 

wherein W is a matrix relating output at source and detector 
positions r s and r d at time t to input at a voxel’s position r, 
A is a regularization matrix, chosen for convenience to be 
diagonal but selected in a way related to the ratio of the 
noise, <nn> to fluctuations in the absorption (or diffusion) X- 
that we are trying to determine: 

A—kjbij with A.y=<nn>/<AXjAXj > 

Y is the data collected at the detectors, and X* is the kth 
iterate toward the desired absorption information. 

A distinct advantage of the present method, over existing 
45 inverse reconstruction methods for imaging objects located 
in highly scattering turbid media using diffusive photons, is 
that the present method can be done with fewer iterations, 
thereby making possible real-time imaging. 

An inverse algorithm, which combines a two dimensional 
50 (2D) matrix inversion with a one -dimensional (ID) Fourier 
transform inversion, is also disclosed herein to obtain 
images of three-dimensional (3D) objects hidden in turbid 
scattering media. To obtain such 3D images, the source - 
detector pairs are arranged in a 2D plane (the x-y plane), 
55 which are scanned along the z direction. Data is obtained 
from S-D pairs at z positions for processing. The reference 
system is assumed to be uniform and infinite along the z 
direction. In the above experimental setting, the weight 
function, W, satisfies translation invariance. In other words, 
60 it is a function of z-z 0 , where z 0 is the z-position of 
source -detector plane and z is z-coordinate of a voxel. 

After making a one -dimensional Fourier transform of W 
over z-z 0 , we obtain K independent 2D matrices, W 2D (k), 
parameterized by k, with K the number of grid-points in the 
65 Fourier k-space. We separately calculate the K inverse 
matrices [W 2D ) 7 W 2Z> (k)+A(k)] _1 , k=l,2,3, . . . ,K, where 
A(k) is a matrix for regularization. These inverse matrices 
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are stored as a database for later reconstruction of the image 
in different hidden object cases. 

Experimental data in time-resolved (or frequency 
domain) are Y(r^, r s , t(oo), z 0 ), where r d and r s are (x,y) 
coordinates of detector and source, t is the time slice (or to 
is frequency). We make a Fourier transform over z 0 to obtain 
Y(r^, r s , t(oo), k), k=l,2,3, . . . ,K. We then calculate the image 
in the k space using the following matrix multiplication: 

Xik^YikyW^kjW^ky W^+A^)]- 1 , k=l,2,3, . . . ,K 

Then, the inverse ID Fourier transform of X(k) over k 
generates the 3D distribution of the change of absorption 
coefficients or scattering coefficients. 

The aforementioned method of 3D imaging using matrix 
inversion and Fourier transform greatly reduces the compu- 
tational burden compared to standard 3D matrix inversion 
methods, while maintaining mathematical integrity in the 
use of Fourier Transform. 

To demonstrate the capacity of the above -described 
method, we used said method to obtain a 3D image with 
28x28x32 voxels. An image was reconstructed by running 
about 30 minutes on a Silicon Graphic Indy 4400 computer 
for 3D. In this manner, the data may be taken and the image 
calculated while a patient waits for results. 

The present invention is also directed to systems con- 
structed in accordance with the foregoing methods. 

In the description which follows, reference is made to the 
accompanying drawings which form a part thereof and in 
which is shown by way of illustration specific embodiments 
for practicing the invention. These embodiments will be 
described in sufficient detail to enable those skilled in the art 
to practice the invention, and it is to be understood that other 
embodiments may be utilized and that structural changes 
may be made without departing from the scope of the 
invention. The following detailed description is, therefore, 
not to be taken in a limiting sense, and the scope of the 
present invention is best defined by the appended claims. 

BRIEF DESCRIPTION OF THE DRAWINGS 

The accompanying drawings, which are hereby incorpo- 
rated into and constitute a part of this specification, illustrate 
various embodiments of the invention and, together with the 
description, serve to explain the principles of the invention. 
In the drawings wherein like reference numerals represent 
like parts: 

FIGS. 1(a) and 1(b) are simplified schematic perspective 
and section views, respectively, of an experimental setup 
used to test the imaging method of the present invention; 

FIG. 2 is a graphic representation of various time- 
resolved curves obtained using the experimental setup of 
FIGS. 1(a) and 1(b), curves (a) through (c) being obtained 
by positioning an object in the three locations (a), (b) and 
(c), respectively, of FIG. 1(b), and the curve labelled “Uni- 
form” being obtained by not positioning the object within 
the sample medium; 

FIGS. 3(a) through 3(e) are comparative images of the 
object of FIG. 1(a) obtained using the method of the present 
invention with different regularization parameters; 

FIGS. 4(a) through 4(A) are comparative images of the 
object of FIG. 1(a) obtained using the method of the present 
invention with different numbers of iterations; 

FIGS. 5(a) and 5(b) are comparative images of the object 
of FIG. 1(a) obtained using the method of the present 
invention with (a) a comparatively great degree of regular- 
ization and no iteration; and (b) a comparatively low degree 
of regularization and 2000 iterations, respectively; 
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FIG. 6 is an image of the object of FIG. 1(a) obtained 
using the conventional Simultaneous Algebraic Reconstruc- 
tion Technique (SART) after 3 hours and 100 iterations; 

FIGS. 7(a) and 7(b) are schematic views of an experi- 
5 mental setup used to generate a 3D image of two objects 
hidden in a scattering medium; 

FIGS. 8(a) and 8(A) are schematic representations of two 
different arrangements of a pair of 3D objects hidden in a 
scattering medium; 

FIG. 9(a) represents a series of images of the objects of 
FIG. 8(a) shown at different places in the z direction; 

FIG. 9(A) represents a series of images of the objects of 
FIG. 8(A) shown at different places in the z direction; 

15 FIG. 10(a) is a simplified schematic diagram of an experi- 
mental arrangement used to implement the imaging method 
of the present invention, the arrangement using an electronic 
time-gating technique with a gated image intensifier and a 
CCD camera; 

20 FIG. 10(A) is a simplified schematic diagram of a sample- 
cell arrangement used to accompany the imaging arrange- 
ment of FIG. 10(a). Also presented is a schematic diagram 
of a sequence of shadowgram-like images that are two- 
dimensional intensity distributions as a function of gate 
25 delay, I(x,y,t). Integrated intensity over the same small area 
of all the shadow images in the sequence plotted as a 
function of gate position generates a temporal intensity 
profile, also displayed in FIG. 10(A) above the sequence of 
image frames. Temporal intensity profiles may also be 
30 obtained by plotting the average intensity over a small area 
or the intensity value at a particular point on the shadow 
images as a function of gate position; and 

FIG. 11 shows a sequence of four shadow images mea- 
sured with the electronic time-gated imaging camera system 
35 of the present invention, the frames (a), (b), (c) and (d) 
corresponding to images taken at gate positions of -50, 0, 50 
and 100 ps, respectively, with respect to the ballistic pulse. 

DETAILED DESCRIPTION OF PREFERRED 
40 EMBODIMENTS 

As indicated above, the present invention is based upon 
the present inventors’ discovery that objects hidden in 
highly scattering turbid media can be imaged by a novel, fast 
45 and noise -resistant inverse method which comprises, among 
other things, transilluminating a highly scattering turbid 
medium using a plurality of intersecting source-detector 
pairs and using intensity measurements of the diffusive 
component of the detected light signals taken over a plural- 
5Q ity of time windows (or “time-slices”) in a novel inverse 
reconstruction algorithm. 

The following discussion provides the theoretical basis 
for the present invention: A scattered signal is, in general, a 
function of key parameters characterizing the medium. The 
55 forward problem is to build a mathematical model that 
allows computation of signals when a spatial distribution of 
key parameters is given. This mathematical model is clearly 
based on the fundamental physical principles governing the 
interaction. The inverse problem is to find the distribution of 
6Q these key parameters that best fit the given signal data. For 
photon propagation in turbid media, the key parameters 
characterizing the interaction are the absorption length, the 
transport mean free path, the scattering mean free path, and 
the index of refraction of media. The parameters are in 

65 general functions of position ( r ) and wavelength (k). In the 
case of optical diffusion tomography, the most commonly 
used key parameters are absorption and transport scattering 
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coefficients, which are the inverse of absorption length and 
transport scattering mean free path, respectively. The math- 
ematical model governing the photon propagation is the 
photonic transport theory or its approximations, such as the 
diffusion theory or the telegrapher’s equation. 

To show the basic principle of this technique, diffusion 
theory is applied here to model light propagation in turbid 
media. The diffusion equation in turbid media is 

[ ^ + Cfi a (r) - V (<?Z)(r)V)j l(r, t) = S(r, t) 


photons that migrate from this voxel to a detector. AV is 
volume of a voxel. Therefore, W in Eq.(3.1) represents a 
normalized effect that photons originating from the source 
migrate through voxel j, then further migrate to the detector 
5 during a total time of t. In the case of diffusion or scattering 
change, both G° in the integral of Eq.(3.1) are replaced by 
VG°. 
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W[(r d , r s , r), rj] = 


(3.2) 


AV 


av n 

: ( drVG 0 (r d , rj, t - r) • V G° (rj, r s , r) 

r d , r s , t)J o 


where I( r ,t) is the field of photon density, S( r ,t) is the 
source of laser light, ^ i a is the absorption coefficient, D=l/ 
[3(a fl +(l-g)a s )] is the diffusion coefficient, with /u s the 
scattering coefficient, g the scattering anisotropic factor and 
c the speed of light in the medium. 

We are interested in the following problem: given inten- 
sity measurements taken around the media either with 

objects located therein, I( r ,t), or without objects located 

therein, I 0 ( r ,t), reconstruct a spatial map of the change in 
optical parameters (absorption and diffusion coefficients) 
due to the presence of hidden objects. 

The forward problem of scattering can be represented in 
a matrix form written as 

Y=WX+n (2) 

where Y has M=(number of S-D pairs)x(number of time 
slices) elements, which is the change in intensity profile and 

is a function of the position of source r the position of 

detector r d , and the time-slice in the profile t. X has 
N=(number of voxel) elements wherein a voxel is a space 

grid element. X y =cA// fl ( r y ) or cAD( r y ) is the (time- 
independent) change in the absorption coefficient or the 
diffusion coefficient in the voxel j due to the presence of 
hidden objects, and n is signal detector noise. W is an MxN 
matrix. Its element is the weight function or “photon mea- 
surement density functions.” This weight function repre- 


15 One preferred geometry for experiments with the elec- 
tronic time-gated imaging camera (TGIC) system, shown in 
FIG. 10(b), requires the use of cylindrical co-ordinate sys- 
tem. The source is located on the origin. The G°(r,r s ) in Eqs. 
(3.1) and (3.2) is c|) independent, and G°(r,r^) depends on 
2Q (<|)-<|) rf ), if a cylinder symmetry and a uniform background is 
assumed. In above experimental setting, the weight function, 
W, satisfies translational invariance. In other words, it is a 
function of (4>-c|>^) 3 where c|)^ is the polar angle of detectors 
and § is the polar angle of a voxel. Therefore, W 3Z> has 
25 indices of p d , t; p, z, (|>-<|)^ while z d is fixed. 

After making a one -dimensional Fourier transform of 
W 3Z? over (<b-<b d ), we obtain K independent 2D matrices, 
parameterized by k, with K the number of grid-points in the 
Fourier k-space. W ^ is function of (p^,t;p,z). Information 
30 on time slices provides an equivalent ID data for inverse of 
image on (p,z) voxels in k space. 

We separately calculate the K inverse matrices [W 2Z> (k) 
7 W 2Z? (k)+A(k)] _1 , k= 1,2,3, . . . ,K, where A(k) is a matrix for 
regularization. These inverse matrices are stored as database 
35 for later reconstruction of image for different hidden objects. 

Experimental time-resolved data are Y(p^,c|)^,t). We make 
a Fourier transform over § d to obtain Y(p^,t,k), k=l,2, 
3, . . . ,K. We then calculate the image in the k space using 
the following matrix multiplication: 

40 

X(k)=Y(k) T W 2D (k)[W 2D (k) T W^(k)+ AC*)- 1 ! k=l,2,3, . . . ,K (3.3) 


sents the change in signal at the detector placed at r d at time 

t from a source signal originating at r 5 associated with a 
unit change in absorption coefficient (or diffusion 

coefficient) at r ■ of the reference medium, as implemented 
by a diffusion process that carries the signal from source at 

r to the detector at r d . Its value is theoretically calculated 
according to the mathematical model, for example, the 
diffusion theory in the case shown below. 

Using Green’s function formula, for the case of absorp- 
tion change, cA ju a , the element of W is given by 


W [(?</, r s , t), rj] = — r drG°(r d , rj,t- r)G°(r y , r s , r) 

G°(r d ,r s ,t)J o 


(3.1) 


where G°( r y , r ( .,t) is the Green’s function of Eq.(l) for a 
reference system (the system with no hidden object) which 

gives the normalized photon intensity at r ■ at time t as 

migrating photons originate from r i at time t=0. In Eq.(3.1), 

the second Green’s function in the integral, G°( r y , r s ,x ) 
represents the photons that migrate from a source to a given 

voxel element. The first one, G°( r d , r ., t-x), represents the 


Then, the inverse ID Fourier transform of X(k) over k 
generates the 3D distribution of the change of absorption or 
45 scattering coefficients on (p,c|),z) voxels. 

According to certain previously -used inversion algo- 
rithms employed in image reconstruction, such as Algebraic 
Reconstruction Technique (ART), Simultaneous Algebraic 
Reconstruction Technique (SART) or Conjugate Gradients 
50 (CGD), the fractional signal change Y=-(I-I 0 )/I 0 , which 
corresponds to the first order perturbation, is used. Such 
inversion reconstruction algorithms, however, generally 
requires a great number of iterations in order for a clear 
image map of hidden objects to be obtained. Such iterations 
55 frequently can take an amount of calculation time that 
exceeds the practical time limit for clinical applications. 

A more accurate non-perturbative approximation based 
on a cumulant expansion has been developed by us, which 
fits the experimental data by Y=-ln(I/I 0 ). This treatment, to 
60 some extent, automatically includes higher order non-linear 
contributions. The inversion is a reconstruction of optical 
parameters from signals. In order to make a fast inversion 
for clinical applications, we have developed the following 
novel inverse reconstruction procedure based on a modifi- 
es cation of the iteration formula described in Shaw, “Improve- 
ment of the Numerical Resolution of An Instrument by 
Numerical Solution of the Integral Equation,” J. Math. 
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Analys. and Applic., 37:83-112 (1972), which is incorpo- 
rated herein by reference: 

X {k+l)T = [Y t W + X {k)T A][W t W + A]" 1 (4) 

where Y, (-ln(I/I 0 )), is the input data for inversion. and 
X(* +1 ) are k th and (k+l) r/l iteration solutions, respectively, for 
the change of parameters. The initial value of X is set to zero 
in the absence of prior information about hidden objects. T 
represents the transpose matrix. The matrix A, which regu- 
larizes or converts an ill-posed problem to a well-posed 
problem, is chosen to be diagonal, A^—X-8^-. X,=<nn>/ 
<AXAX>, with < . . . > the statistical average, and AXy= 
X--<X>. Thus, Xj represents the ratio of mean square of 
noise to the expected mean square deviation in the change of 
parameters (structural noise). 

Since W does not relate to the position and optical 
parameters of hidden objects, having W for a given reference 
system (either homogeneous or inhomogeneous), the inverse 
matrix [W^W+A] -1 can be precalculated and stored as a 
database, which can be used in clinical cases to find hidden 
objects. This strategy separates time-consuming work in 
inversion processes into two parts. Pre-computation is time- 
consuming and may require the use of a supercomputer. 
However, it greatly reduces the burden of performing 
inverse processes in daily clinical cases. 

One of the main difficulties in inversion is that most 
inverse problems are ill-posed. If no noise existed, an 
inverse image could be reconstructed from Eq.(2) via: 
X r =Y r W[W r W ] _1 . The existence of a small amount of 
noise should not have much effect on the results obtained. 
However, the matrix W r W is generally ill -posed or near- 
singular. In other words, this matrix has some near-zero 
eigen values. This occurs when the number of measurements 
M is less than the number of unknown parameters N or when 
the measured signals are not fully independent of each other. 
In this case, the effect of small noise will be magnified in the 
inverse process, and the solution of inversion becomes 
unstable. One of the methods to overcome the ill -posed 
problems is regularization, which adds a A to the matrix 
W r W. A is relatively small compared to the matrix W r W; 
therefore, regularization makes the matrix well-posed or 
well-defined, but does not change the characteristics of the 
image map. The value of A is set according to the experi- 
mental noise level and structural noise present in the change 
of optical parameters. If ||AW r W||<l, one may use the 
following expansion to correct the error introduced by 
regularization: 

[W r W\~ 1 =[B-A\~ 1 =B~ 1 +B~ 1 AB~ 1 +B~ 1 AB~ 1 AB~ 1 + ... (5) 

The iterations in Eq.(4) are equivalent to the above expan- 
sion. 

If noise n exists, Y=WX+n. One wants to find an inverse 
operator, X=LY, which gives the estimated value X of X, 
that, in the sense of statistical average, has minimum error. 
Shaw has derived the following solution: 

Jt=[Y T W<AXAX>+<X><nn>]x[<AXAX>W T W+<nn>]~ 1 (6) 

with < . . . > statistical average, AX=X-<X>. We see that 
Eq.(6) is a well-posed expression. This approach is the 
so-called “well-posed stochastic extension of ill-posed lin- 
ear processes” or “statistical regularization.” Comparing 
Eq.(4) and Eq.(6), A in Eq.(4) is related to <nn>/<AXAX>. 
Thus, Xj represents the ratio of mean square of detector noise 
to the expected mean square deviation in the change of 
parameters (structural noise). 
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The key point is that since W is only reference -related, 
and does not relate to the position and optical parameters of 
hidden objects, the inverse matrix [W^W+A ] -1 of a given 
pattern of the reference system, which can be either homog- 
5 enous or inhomogeneous, can be pre-calculated and pre- 
stored as a database for various patients for clinical use. 

Our tests show that this novel inverse procedure con- 
verges after only 2-4 iterations, which is much faster than 
the ART or SART procedures currently used. It takes about 
10 one minute to form an image in our reconstruction tests on 
a Silicon Graphic Instruments (SGI) workstation computer. 

To demonstrate the present technique, the following 
experiment was performed. Schematic drawings of the 
experimental arrangement are shown in FIGS. 1(a) and 1(b). 
15 Ultrashort pulses of 100 fs pulse duration at a wavelength of 
625 nm were coupled into a rectangular shape random 
medium by an optical fiber. This geometry was chosen 
because it is easy to construct an analytical expression of 
Green’s function G° in a finite sized 3D volume, which is 
20 suitable for testing the physical modeling and the inversion 
algorithm. The scattered pulses around the medium were 
collected by 7 optical fibers and time resolved by a streak 
camera. The size of the sample was 60x60 mm and the 
height was 90 mm. The fibers were located at mid-plane. 
25 The source fiber was at the center of the boundary wall. Two 
fibers each were placed on the opposite and the side bound- 
ary walls of the sample, and one fiber detector was placed on 
the incident wall. The scattered pulses collected by the 7 
detection fibers and a reference pulse were coupled onto the 
30 input slit of the streak camera for multi-channel detection. 
The medium was Intralipid solution with a transport mean 
free path of 2.5 mm and absorption length of about 500 mm. 
The hidden object was a black painted aluminum rod of 
diameter of 4.8 mm. Temporal profiles of the scattered 
35 pulses around the medium were measured. In this experi- 
mental setting, the structure was z-independent, so the 
reconstructed image is presented as a 2-D image. However, 
the present method could also be used as a 3-D imaging 
model. 

40 As an example, a typical set of temporal profiles measured 
for one S-D pair is shown in FIG. 2. In this case, the detector 
is located straight across the medium in the line of incidence. 
Measurements of the uniform medium without the object 
and with the object placed in three locations 90 degrees apart 
45 are also shown. It is clearly shown that the intensity is 
reduced by introducing the absorbing object, and a larger 
reduction was observed when the object was placed in the 
line of the source detector pair. The measurement of 3 object 
positions is equivalent to 3 measurements by rotating the 
5 q S-D pairs by 90 degrees while fixing the object location. 

In the image reconstruction, an area in the sample was 
divided into 28x28=784 voxels, which corresponds to an 
actual voxel size of 2x2 mm. Measurements of 7 S-D pairs 
with the above mentioned object locations were used. For 
55 each temporal intensity profile, intensity data at 40 time 
slices uniformly distributed from 602 ps to 2045 ps were 
used. Thus, the total number of input data points was 
21x40=840. This is one of the important features of our 
method — the use of hundreds of time slices from the tem- 
60 poral profile of one S-D pair for image reconstruction gives 
much more information than data obtained in the frequency 
domain with much higher speed resolution in sub -mm to 
mm range to detect a defect. The robustness of the method 
was demonstrated with a small number of fibers. In our 
65 programs, Y and W are magnified by a factor of 100. 

Our regularization parameters are spacing independent: 

A (rj)=ko exp (dk\r r r cent J) (7) 
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FIGS. 3(a) through 3(e) show images obtained for in 
Eq.(4) with different regularization parameters. Using large 
X 0 , the image of the hidden object has already appeared after 
the first iteration. With a decreasing of X 0 down to 10~ 2 , the 
image is gradually broken. 

FIGS. 4(a) through 4(h) show the results of different 
numbers of iterations, when the small regularization X 0 is 
used. We see that after m=50 iterations, the image of the 
hidden object begins to appear. After m=500, the image 
becomes stable as no visible change in the shape of the 
object occurs between m=500 to m=2000. 

FIGS. 5(a) and 5(b) show an image obtained where (a) 
X o =1000, without iteration, by running 1 minute on a Silicon 
Graphic Incorp. (SGI) Indy computer where [W^W+A] -1 is 
stored; and (b) X o =10 -2 , with 2000 iterations, by running 3 
hours on the same computer, respectively. As can be seen, 
the two images are similar, but the latter has better resolu- 
tion. Noting that the difference of regularization parameters 
for these two images is up to 5 orders of magnitude, the 
stability of the imaging process for the present technique is 
confirmed. 

In both FIGS. 5(a) and 5(b), the presence of an absorber 
can be clearly seen. The location of the maximum change is 
about 7 mm away from the center which is about the actual 
location of the object (10 mm). The full width of half 
maximum of the absorber is 8 mm which is comparable to 
the diameter of the object (4.8 mm). The fluctuation at other 
spatial locations are less than one third of the peak value. 
Using different regularization parameters, the same location 
of the absorber was found, but the maximum peak value and 
the full width of half maximum were different. One possible 
reason that the center of the image deviates from that of the 
experimental set is that the hidden object is not a pure 
transparent scattering absorber, and the possible effect of 
reflection on the surface of the rod was not included in our 
physical model. 

For comparison, the result using the SART inversion 
algorithm, commonly used by others, after 100 iterations 
(taking X=0 as initial input), which runs 3 hours on the same 
SGI workstation, is shown in FIG. 6. There is no clear image 
of the absorber, noticing that the maximum absorption 
coefficient is only 0.0025 which is about 40 times smaller 
than that obtained by our new algorithm. 

We have also tested our technique on simulated temporal 
data with 15% white noise added. A clear image of a single 
absorber on the correct position was obtained. These results 
show that our approach is noise-resistant. This inversion 
method can include a change of diffusion coefficients, AD( 

r ) or A ju s ( r ), as described before for detecting scattering 
defects in a highly scattering host medium. 

It should be noted that a reconstruction algorithm which 
is suitable for clinical applications must be reasonably rapid 
in running time while maintaining mathematical integrity. 
An inverse algorithm mainly based on Fourier transform 
runs in the fastest way. In order to make a Fourier transform 
deconvolution available, the condition of “translation invari- 
ance” should be satisfied. In other words, the propagator 
K(y,x) must be the function of y-x. This condition is violated 
in experimental settings in the following three ways: (1) the 
media is not uniform; (2) there are boundaries on which the 
light sources and detectors are setting; (3) the propagator 
here is a function of (r s ,r^,r); therefore, it is impossible to 
make a Fourier transform for both r s -r and r^-r. Among 
above listed violations, (1) is relatively soft, since we can 
assume a reference with uniform background, and put non- 
uniform as an image. On the other hand, reconstruction 
based on the inverse matrix method does not require the 
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condition of translation invariance and can be used for more 
general cases. For purposes of clinical applications, a three- 
dimensional (3D) reconstruction algorithm is necessary. 
Because light diffuses in 3D space, the measurements on a 
5 specific layer are affected by the hidden objects on other 
layers. In extending our approach to the 3D case, we face a 
new difficulty. The computational complexity to inverse a 
matrix is about N 3 , where N is the size of the matrix. If one 
supposes that the size on the third dimension is divided into 
10 30 grids, the computing complexity to inverse a matrix in the 
3D case is about 27,000 times that in the 2D case, which 
makes the precomputation of an inverse matrix for 3D case 
unfeasible. To overcome this difficulty, we propose adding a 
ID Fourier transform to the 2D matrix inversion. 

15 Therefore, a novel inverse algorithm for obtaining images 
of 3D objects hidden in scattering media is taught in this 
invention, the method combining a 2D matrix inversion with 
a ID Fourier transform inversion. In our approach the 
conditions for Fourier transform available are properly 
20 treated. This method greatly reduces the computational 
burden, compared to standard 3D matrix inversion methods. 
This approach can be used for both time-resolved data and 
frequency domain data. Using the above mentioned pre- 
computation and database storage strategy, we tested the 
25 present technique. It takes about 30 minutes on a Silicon 
Graphic Indy 4400 computer to obtain a 3D image of 
28x28x32 voxels of mm objects in cm’s scattering media. 

The source -detector pairs are arranged on the boundary of 
scattering media, such as breast and brain, in a 2D plane (the 
30 x-y plane), which is scanned along the z direction. 
Accordingly, z coordinates of sources and detectors are the 
same, thereby eliminating the violation (3) mentioned 
above. The reference system is assumed to be uniform and 
infinite along the z direction. This assumption leads to some 
35 deviation in the reconstructed image due to finite z-size in 
experimental setting. Since there is no source and detector 
setting on the z-end plane, this end effect becomes less 
serious. Under the above experimental setting and 
assumption, W satisfies z-translation invariance. In other 
40 words, it is a function of z-z 0 , where z 0 is z-position of 
source -detector plane and z is z-coordinate of a voxel. After 
making a one -dimensional Fourier transform over z-z 0 we 
obtain K independent 2D matrices, W 2D (k), parameterized 
by k, with K the number of grid-points in the Fourier 
45 k-space. We separately calculate the K inverse matrices 
[ W 2 o(k) 7 W 2jD (k) + A(k)]“ 1 , k=l,2, . . . ,K, where A(k) is a 
matrix for regularization. These inverse matrices are stored 
as database for later reconstruction of image in different 
hidden object cases. The computational complexity now is K 
50 times that of inverse of a W^ matrix, which is much less 
than that of inverse of a W 3Z> matrix. 

The experimental data in time-resolved (or frequency 
domain) are Y(r^,r s ,t(oo),z 0 ), where r d and r s are (x,y) coor- 
dinates of detector and source, respectively, and t is the slice 
55 (or co is frequency). We make a Fourier transform over z 0 to 
obtain Y(r^,r s ,t(co),k), k=l,2,3, . . . ,K. We then obtain the 
image in the k space using the following matrix multiplica- 
tion: 

JST(A:)=n*) r W2^(*)[ 7 'W 2 ^(A;)+ A(*)]- \ k=l,2,3, . . . ,K 
6U 

Then, the inverse ID Fourier transform of X(k) over k 
generates the 3D distribution of the change of absorption 
coefficients or scattering coefficients. 

A simulated experimental setting for a 3D image obtained 
65 using a plurality of source -detector (S-D) pairs is shown in 
FIGS. 7(a) and 7(b). The S-D pairs are arranged around a 
square of size 60x60 mm 2 . This arrangement of S-D pairs is 
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then rotated 90 degrees and 180 degrees and the plane is 
scanned 32 steps along z in a range of 120 mm. For each 
simulated temporal profile, intensities at 40 time slices 
uniformly distributed from 602 ps to 2045 ps are taken. The 
number of sampling data are 7x3x40x32. Data both with 5 
and without hidden objects are computed to produce 
Y=-1001n(I/I o ), adding uniform distributed noise in ±10% 
range, where I and I 0 are intensity with and without hidden 
objects, respectively. The volume of sample is divided into 
28x28x32 voxels. Two arrangements of hidden absorption 1Q 
objects are separately shown in FIGS. 8(a) and 8(b). The 
cross-sections of the hidden objects are 4x4 mm 2 and 6x6 
mm 2 . In case (a), both objects are located from the 9th 
z-layer to the 13th z-layer (about 20 mm in length). In case 
(b), one is located from the 9th layer to the 13th layer (about 15 
20 mm in length) and another is located from 17th layer to 
the 20th layer (about 16 mm in length). The W is calculated 
using Green’s function method based on the diffusion equa- 
tion. For absorbing coefficients, it is given by the following 
equation: 2Q 


W 3 D(r s ,r d , r,z-Zo,t) = 


(4 nDct) i!2 
G%(r s / d , t) 


£ 


drG\ D {r d , r, t-r)G\ D {r, r s , r)x 


[4 nDc(t-T)\ 1/2 [4 ttDct] 1/2 exp| 


t(z-Zo ) 2 ' 
ADc{t - t)t ^ 


25 


Where G^ 0 is 2D Green’s function for reference, D is the 30 
diffusion constant and c is the speed of light in the medium. 
Scattering objects with different \ t can be also calculated as 
mentioned above for 2D case. The background transport 
mean free path is chosen as l r °=2.5 mm, and the absorption 
length l a °=500 mm. The absorbing objects have l fl =2.5 mm. 35 
In making Fourier transform and computing inverse matrix 
[W 20 (k) :r W 2D (k)+A(k)] _:L , the regularization matrix is set as 
A(k)=(X+pk 2 )I, with \ and (3 are parameters for regulariza- 
tion. In this calculation X=10 2 and p=0 are taken. The images 
of the absorption objects are separately shown in FIG. 9(a) 40 
and FIG. 9(b), which correspond to FIGS. 8(a) and FIG. 
8(b), respectively. In FIGS. 9(a) and 9(b), index K is the 
index of layer in z coordinate. The images of hidden objects 
appear in the right layers where the objects are arranged, 
with a resolution of about 4 mm. There are some aliased 45 
images appearing near the z-end because of discrete sam- 
pling. 

The experimental arrangement for implementing the 
present technique when using electronic time-gated imaging 
camera (TGIC) system is shown schematically in FIGS. 50 
10(a) and 10(b). Ultrashort light pulses of approximately 
100 fs duration from a Ti:sapphire laser and amplifier system 
operating at a wavelength of 800 nm and a repetition rate of 
1 kHz are launched into the scattering medium in which the 
object is embedded. One preferred shape of the random 55 
medium is a cylindrical disc. The radius of the disc is 
preferably much larger than its thickness. The cylindrical 
geometry is chosen because it is easier to construct an 
analytical expression for Green’s function for testing the 
physical model and inversion algorithm. 60 

The emergent pulse from the other end face of the 
cylindrical sample holder is amplified and different temporal 
slices of the pulse are recorded by the electronic time -gated 
image intensifier and CCD camera system. The output of the 
TGIC system is a shadowgram-like image, or more 65 
generally, a two-dimensional time-dependent intensity dis- 
tribution I(x,y,t). These are shown schematically in FIG. 


10(b) as sequence of frames. Integrated intensity over the 
same small area of all the shadow images in the sequence 
plotted as a function of gate position generates a temporal 
intensity profile, also displayed in FIG. 10(b) above the 
sequence of image frames. Temporal intensity profiles may 
also be obtained by plotting the average intensity over a 
small area or the intensity value at a particular point on the 
shadow images as a function of gate position. I(x,y,t) may be 
readily converted to Y(p^,c|)^,t) for use in the inversion 
algorithm. 

The sample we used was a glass slide on which scattering 
particles of size varying from 0.1-1 .5 mm were attached. It 
was placed at the center of the 5.5 cm thick sample cell. The 
scattering medium was a 1% solution of Intralipid 10% in 
water. A sequence of four frames taken at gate positions of 
-50, 0, +50 and +100 ps are displayed in FIG. 11 as frames 
(a), (b), (c) and (d), respectively. The zero time corresponds 
to the arrival time of the ballistic light that is determined by 
performing the experiment using clear water in the sample 
cell. The sequence of frames clearly demonstrates that at 
later times the diffusive photons blur out the shadow image, 
as in frame (d). For highly-scattering thick samples, such as 
human breast, blurred images as in frame (d) or worse are 
expected. The present invention provides an algorithm to 
reconstruct images from such blurred shadows. The 
sequence of temporal slices, that is, the shadowgrams taken 
at different time intervals enables determination of location 
of the object inside the scattering medium. 

Scattering light pulse from the side wall of the cylindrical 
cell may be collected as well, as shown in FIG. 10(b), to 
provide additional intensity profiles for image reconstruc- 
tion. Scattered light pulse from the side wall may be 
recorded using the same camera system, or preferably 
another identical camera system or systems. Sample cells of 
other geometries, such as a cube, a rectangular 
parallelipiped, or of an arbitrary shape may also be used with 
corresponding modifications of image reconstruction algo- 
rithm. 

The following comments, observations, objects, features, 
uses, applications and/or advantages may be made about the 
present invention: 

(1) Temporal intensities extracted at multiple time slices 
of scattered light profile measured using multiple 
source detector pairs around turbid media are used as 
input data for the inversion to obtain the image maps of 
the media. 

(2) An image reconstruction algorithm for imaging highly 
scattering turbid media is developed. The inversion 
procedure consists of input scattered light intensity data 
at different time slices, physical modeling of light 
propagation in random media, inversion computation 
algorithm, reconstructed image of spatial distributions 
of key optical parameters of random media, and image 
displays. 

(3) Using different theoretical models for describing pho- 
ton migration in random media for imaging. The theory 
includes transport theory for photon migration, and its 
approximations, such as the diffusion theory and the 
telegraph equations. 

(4) Key optical parameters, such as the absorption 
coefficient, the index of refraction, the transport scat- 
tering length, the scattering length, and the diffusion 
constant, are mapped for imaging the internal structure 
of the turbid media at different light wavelengths. A 
difference in both absorption and scattering can be 
imaged to form a map of the internal structure of the 
turbid media in 3-D. 
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(5) A fast inversion algorithm is developed based on 
Shaw’s principle, Eq. (4), in which the weight function 
matrix, W (Eq.(3)), which relates a change in the 
properties of turbid media to the changes in measured 
light intensities, can be pre-computed and stored for 5 
fast image inversion for each patient. 

( 6 ) A cumulant formula is developed to relate the experi- 

mentally measured intensities to the prediction from 
theoretical models based on properties of the turbid 
media. 10 

(7) 3-D images of spatial distribution of key optical 

parameters of turbid media can be obtained using this 
algorithm. 2-D slices of tomographic images can be 
obtained by taking a section of a 3-D image and 
displaying it to highlight the location of tumors. 15 

( 8 ) Use the algorithm to image non-invasively tumors and 
abnormality in human body, such as tumor growths in 
breasts, tumors in brain, cerebral hemorrhage in brain, 
hemorrhage in internal organs, prostate. 

(9) Using the algorithm with temporal data at various NIR 20 
wavelengths from 700 to 1500 nm for image recon- 
struction for monitoring metabolic states of a given part 

of human body in vivo. 

(10) This novel inverse reconstruction method can be 
combined with shadowgram type of image techniques 25 
to further improve spatial resolution and enhance the 
presence of the abnormality. 

(11) Use single source optical fiber and multiple detector 

fibers for collecting temporal profiles of the scattered 
pulses. 30 

(12) Use multiple input source fibers and multiple detec- 
tor fibers for collecting temporal profiles of the scat- 
tered pulses. A time sharing beam deflector will be 
incorporated to switch the input from one source fiber 35 
to another. 

(13) Multiple S-D optical fibers are mounted on bra-like 
or hat-like holders surrounding a breast or a head 
(brain) to detect temporal intensity profiles of scattered 
pulses from breast or head (brain) for imaging the 4Q 
internal structure of optical properties. 

(14) The lasers used for imaging are Cr 4+ Forsterite, Cr 4+ 
YAG, semiconductor lasers, Nd:YAG, and/or Ti: Sap- 
phire lasers. The laser pulses at different wavelengths 
will be coupled to optical fiber by a time sharing 45 
scheme to deliver to the scattering media for spectral 
imaging. 

(15) Temporal intensity profiles of scattered pulses around 
scattering media are measured using time resolved 
detection methods, including streak camera, ultrafast 50 
pin diodes, fast photomultipliers. 

(16) Time-resolved shadowgram -type images and their 

intensity profiles obtained using the electronic time- 
gated image intensifier and CCD camera system will 
provide a plurality of profiles. For example, the TGIC 55 
system we presently use includes a 384x284 pixel CCD 
camera. Each pixel in this arrangement may be looked 
upon as a detector fiber in the one-source-7-detector 
arrangement that has been described above. Thus, the 
data collection time is shortened and the use of a 60 
plurality of profiles increase the spatial resolution of 
reconstructed image. Since the shadow images are 
taken sequentially in time the depth information and 
hence the location of the object inside the medium is 
obtained. 65 

(17) Reconstructed images are to be displayed in 3-D or 
in 2-D tomography slices. 
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(18) Manipulating images obtained from claim 9 at dif- 
ferent wavelengths to construct difference, sum, ratio 
and/or sum of ratio image maps at two wavelengths or 
several wavelengths to highlight tumors and to differ- 
entiate benign from malignant tumors. 

The embodiments of the present invention described 
above are intended to be merely exemplary and those skilled 
in the art shall be able to make numerous variations and 
modifications to it without departing from the spirit of the 
present invention. All such variations and modifications are 
intended to be within the scope of the present invention as 
defined in the appended claims. 

What is claimed is: 

1. A method of imaging an object located in a highly 
scattering turbid medium, the high scattering turbid medium 
being an in vivo biological tissue, said method comprising 
the steps of: 

(a) illuminating the object through the highly scattering 
turbid medium with a pulse of light, the light emergent 
from the highly scattering turbid medium consisting of 
a ballistic component, a snake-like component and a 
diffusive component; 

(b) determining the intensity of said diffusive component 
at a plurality of points in time; and 

(c) using said intensity determinations to form an image 
of the object in the highly scattering turbid medium, 
said using step comprising using a mathematical inver- 
sion algorithm to generate an image of the highly 
scattering turbid medium, said mathematical algorithm 
being 

X {k+l)T = [ y t w + x {k)T a][w t w + A ] -1 

wherein W is a matrix relating output at source and 
detector positions r s and r d , respectively, at time t to 
position r, A is a regulation matrix, chosen for conve- 
nience to be diagonal but selected in a way related to 
the ratio of the noise, <nn> to fluctuations in the 
absorption (or diffusion) X- that we are trying to 
determine: 

A—Xjbij with Ay=<nn>/<AXjAXj> 

Y is the data collected at the detectors, and X* is the kth 
iterate toward the desired absorption information. 

2. The method as claimed in claim 1 wherein said pulse 
of light has a wavelength in the visible to near infrared 
region of the spectrum. 

3. The method as claimed in claim 2 wherein said pulse 
of light has a wavelength of about 700 nm to about 1500 nm. 

4. The method as claimed in claim 1 wherein said pulse 
of light is an ultrashort pulse of light. 

5. The method as claimed in claim 1 wherein said pulse 
of light is a laser pulse emitted from a laser selected from the 
group consisting of Ti: Sapphire lasers, Cr 4+ Forsterite lasers, 
Cr 4+ YAG lasers, semiconductor lasers and Nd:YAG lasers. 

6. The method as claimed in claim 1 wherein said in vivo 
biological tissue is selected from the group consisting of in 
vivo human breast tissue, in vivo human brain tissue, in vivo 
human neck tissue and in vivo human prostate tissue and 
wherein said object is a tumor. 

7. The method as claimed in claim 1 further comprising 
the step of detecting cancerous tumors in said in vivo 
biological tissue using fluorescence spectroscopy. 

8 . The method as claimed in claim 1 wherein said 
determining step comprises measuring the intensity of light 
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at a plurality of points in the time frame spanning from 50 
ps to 10 ns after illumination. 

9. The method as claimed in claim 8 wherein said 
determining step comprises measuring the intensity of light 

at hundreds of points in the time frame spanning from 50 ps 5 
to 10 ns after illumination. 

10. The method as claimed in claim 9 wherein said 
hundreds of points are determined by slicing the time frame 
into equal portions. 

11. A method of forming a map of a highly scattering 10 
turbid medium, the high scattering turbid medium being an 

in vivo biological tissue, said method comprising the steps 
of: 

(a) illuminating the highly scattering turbid medium with 

a first pulse of light along a first axis of incidence, 15 
whereby the light emergent from the highly scattering 
turbid medium due to said first pulse of light consists of 
a ballistic component, a snake component and a diffu- 
sive component; 

(b) determining, at a plurality of times, the intensity of 20 
said diffusive component of the first pulse of light 
emergent from the highly scattering turbid medium at 

a plurality of locations; 

(c) illuminating the highly scattering turbid medium with 25 
a second pulse of light along a second axis of incidence, 
said second axis of incidence intersecting with said first 
axis of incidence, whereby the light emergent from the 
highly scattering turbid medium due to said second 
pulse of light consists of ballistic component, a snake 3Q 
component and a diffusive component; 

(d) determining, at a plurality of times, the intensity of 
said diffusive component of the second pulse of light 
emergent from the highly scattering turbid medium at 

a plurality of locations; and 35 

(e) using the intensity determinations from steps (b) and 

(d) to generate a map of the highly scattering turbid 
medium, said using step comprising using a mathemati- 
cal inversion algorithm to generate said map, said 
mathematical algorithm being 40 

X {k+l)T = [Y T W + X {k)T \][W T W + A]~ i 

wherein W is a matrix relating output at source and 45 
detector positions r s and r d , respectively, at time t to 
position r, A is a regulation matrix, chosen for conve- 
nience to be diagonal but selected in a way related to 
the ratio of the noise, <nn> to fluctuations in the 
absorption (or diffusion) X y - that we are trying to 50 
determine: 

A with Xj=<nn>l <AXjAXj > 

Y is the data collected at the detectors, and X* is the kth 55 
iterate toward the desired absorption information. 

12. The method as claimed in claim 11 wherein each of 
said first and said second pulses of light has a wavelength in 
the visible to near infrared region of the spectrum. 

13. The method as claimed in claim 11 wherein each of 60 
said first and said second pulses of light has a wavelength of 
about 700 nm to about 1500 nm. 

14. The method as claimed in claim 11 wherein each of 
said first and said second pulses of light is an ultrashort pulse 

of light. 65 

15. The method as claimed in claim 11 wherein each of 
said first and said second pulses of light is a laser pulse 
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emitted from a laser selected from the group consisting of 
Ti: Sapphire lasers, Cr 4+ Forsterite lasers, Cr 4+ YAG lasers, 
semiconductor lasers and Nd:YAG lasers. 

16. The method as claimed in claim 11 wherein said in 
vivo biological tissue is selected from the group consisting 
of in vivo human breast tissue, in vivo human brain tissue, 
in vivo human neck tissue and in vivo human prostate tissue 
and wherein said object is a tumor. 

17. The method as claimed in claim 11 wherein said 
determining steps each comprise time-resolving light emer- 
gent from the highly scattering turbid medium. 

18. The method as claimed in claim 11 further comprising 
the steps of: 

(f) illuminating the highly scattering turbid medium with 
a third pulse of light along a third axis of incidence, said 
third axis of incidence intersecting a plane defined by 
said first and second axes of incidence, whereby the 
light emergent from the highly scattering turbid 
medium due to said third pulse consists of a ballistic 
component, a snake component and a diffusive com- 
ponent; and 

(g) determining, at a plurality of times, the intensity of 
said diffusive component of the third pulse of light 
emergent from the highly scattering turbid medium at 
a plurality of locations; 

(h) wherein said using step comprises using the intensity 
determinations from steps (b), (d) and (g) to generate a 
map of the highly scattering turbid medium. 

19. The method as claimed in claim 11 further comprising 
the step of detecting cancerous tumors in said in vivo 
biological tissue using fluorescence spectroscopy. 

20. The method as claimed in claim 11 wherein said 
determining step comprises measuring the intensity of light 
at a plurality of points in the time frame spanning from 50 
ps to 10 ns after illumination. 

21. The method as claimed in claim 11 wherein said 
determining step comprises measuring the intensity of light 
at hundreds of points in the time frame spanning from 50 ps 
to 10 ns after illumination. 

22. The method as claimed in claim 21 wherein said 
hundreds of points are determined by slicing the time frame 
at regular intervals. 

23. The method as claimed in claim 11 wherein said map 
is an absorption map. 

24. The method as claimed in claim 11 wherein said map 
is a scattering map. 

25. An apparatus for imaging an object located in a highly 
scattering turbid medium, said apparatus comprising: 

(a) means for illuminating the object through the highly 
scattering turbid medium with a pulse of light, the light 
emergent from the highly scattering turbid medium 
consisting of a ballistic component, a snake-like com- 
ponent and a diffusive component; 

(b) means for determining the intensity of said diffusive 
component at a plurality of points in time; and 

(c) means for using said intensity determinations to form 
an image of the object in the highly scattering turbid 
medium, said using means comprising a computer 
programmed to perform a mathematical inversion 
algorithm, said mathematical algorithm being 

X {k+l)T = [Y t W + X {k)T A ][W t W + A ]" 1 

wherein W is a matrix relating output at source and 
detector positions r s and r d , respectively, at time t to 
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position r, A is a regulation matrix, chosen for conve- 
nience to be diagonal but selected in a way related to 
the ratio of the noise, <nn> to fluctuations in the 
absorption (or diffusion) X- that we are trying to 
determine: 5 

A ij=kfiij with Xj=<nn>l <AXjAXj > 

Y is the data collected at the detectors, and X* is the kth 
iterate toward the desired absorption information. io 

26. The apparatus as claimed in claim 25 wherein said 
illuminating means comprises a plurality of light sources 
and wherein said determining means comprises a plurality of 
light detectors. 

27. The apparatus as claimed in claim 26 wherein 7 to 500 15 
light source/light detector combinations are used to generate 

a 2 to 10 mm map of a breast tissue 5 to 10 cm thick. 

28. A method of forming a map of a highly scattering 
turbid medium, said method comprising the steps of: 

(a) illuminating the highly scattering turbid medium with 20 
a first pulse of light along a first axis of incidence, 
whereby the light emergent from the highly scattering 
turbid medium due to said first pulse of light consists of 

a ballistic component, a snake component and a diffu- 
sive component; 25 

(b) determining, at a plurality of times, the intensity of 
said diffusive component of the first pulse of light 
emergent from the highly scattering turbid medium at 
a plurality of locations; 

(c) illuminating the highly scattering turbid medium with 30 
a second pulse of light along a second axis of incidence, 
said second axis of incidence intersecting with said first 
axis of incidence, whereby the light emergent from the 
highly scattering turbid medium due to said second 
pulse of light consists of ballistic component, a snake 35 
component and a diffusive component; 

(d) determining, at a plurality of times, the intensity of 
said diffusive component of the second pulse of light 
emergent from the highly scattering turbid medium at 4Q 
a plurality of locations; and 

(e) using the intensity determinations from steps (b) and 
(d) to generate a map of the highly scattering turbid 
medium, said using step comprising using a mathemati- 
cal inversion algorithm to generate said map, said 45 
mathematical inversion algorithm combining a two- 
dimensional matrix inversion with a one -dimensional 
Fourier transform inversion. 

29. The method as claimed in claim 28 wherein said first 
and second pulses of light are emitted from respective 50 
sources, wherein said determining step is performed with 
detectors located at said plurality of locations, said source 
and detectors being arranged in a 2D plane (the x-y plane), 
which is scanned along the z direction, the reference system 
being assumed to be uniform and infinite along the z 55 
direction and wherein, under above experimental setting and 
assumption, satisfies translation invariance in z coordinate 
and is a function of z-z 0 , where z 0 and z are the z-positions 

of said 2D plane, and of a voxel, respectively. 

30. The method as claimed in claim 28 further comprising 60 
the steps of: 

(a) making a one -dimensional Fourier transform over z-z 0 

to obtain K independent 2D matrices, W 2Z3 (k), param- 
eterized by k, with K the number of grid-points in the 
Fourier k-space; 65 

(b) calculating the inverse matrices [W 2jD (k) 7 W 2jD (k)+A 
(k)] -1 , k=l,2,3, . . . ,K, where A(k) is a matrix for 
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regularization, said inverse matrices being stored as a 
database for later reconstruction of image in different 
hidden object cases; 

(c) obtaining experimental data in time-resolved domain 
(or in frequency domain), Y(r^,r s ,t(oo),z 0 ), where r d and 
r s are (x,y) coordinates of detector and source, t is time 
slice (or 00 is frequency), and making a Fourier trans- 
form over z 0 to obtain Y(r^,r s ,t(co),k), k=l,2,3, . . . ,K; 

(d) using the following matrix multiplication to obtain K 
images in Fourier k-space: 

X{k)=Y{k) T W 2D {k)[ W^^W^^+A^)]- 1 , k=l,2,3, . . . ,K; and 

(e) making an inverse ID Fourier transform of X(k) over 
k to obtain the 3D distribution of the change of absorp- 
tion coefficients or scattering coefficients. 

31. The method as claimed in claim 30 wherein a multi- 
CPU computer is used to compute in parallel K components 
of X(k) in step (d) of claim 30. 

32. The method as claimed in claim 28 wherein said 
mathematical inverse algorithm is used for time-resolved 
data and frequency domain data. 

33. The method as claimed in claim 28 wherein said 3D 
image of both scattering coefficients and absorption coeffi- 
cients of hidden objects is obtained through said mathemati- 
cal inverse algorithm. 

34. The method as claimed in claim 28 wherein said 3D 
image of scattering coefficients of hidden scattered objects 
in scattering media of different scattering length is obtained 
through said mathematical inverse algorithm. 

35. The method as claimed in claim 29 wherein the shape 
of the highly scattering medium is cylindrical, or rectangular 
or that with other shape of (x-y) cross section. 

36. The method as claimed in claim 29 wherein the 
scanning steps in z direction are 16 to 128 to generate an 
image with 16 to 128 layers in z direction. 

37. The method as claimed in claim 29 wherein the 
sources and detectors are fixed on different layers and 
wherein a computer-control switch is used to shine the 
sources sequentially. 

38. The method as claimed in claim 29 wherein the 
sources and detectors are located on the boundary of the 
highly scattering turbid medium or inside the highly scat- 
tering turbid medium. 

39. The method as claimed in claim 29 wherein the highly 
scattering turbid medium is an in vivo biological organ 
selected from the group consisting of human breast, human 
brain, human neck, human prostate, heart, limb, cervix and 
kidney. 

40. The method as claimed in claim 1 wherein the detector 
is an area camera, such as a CCD, and equipped with means 
to record temporal slices of the scattered pulse, such as, a 
gated image intensifies 

41. The method as claimed in claim 40 wherein the gated 
image intensifier has a gate width of 80 ps and gate position 
can be changed in steps of 25 ps. 

42. The method as claimed in claim 40 wherein the in vivo 
biological tissue is from a human breast, said in vivo 
biological tissue being disposed in a cell sample holder, and 
wherein said cell sample holder contains in addition to the 
in vivo biological tissue an index-matched fluid of similar 
scattering characteristics. 

43. The method as claimed in claim 42 wherein the cell 
sample holder is cylindrical. 

44. The method as claimed 42 wherein the cell sample 
holder has a cube-like, sphere-like or any other geometrical 
shape. 
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45. The method as claimed in claim 40 wherein the in vivo 
biological tissue is disposed in a cell sample holder and 
wherein the scattered light is collected from end faces and/or 
side faces of the cell sample holder. 

46. The method as claimed in claim 45 wherein the gated 
image intensifier has a gate width of 80 ps and gate position 
can be changed in steps of 25 ps. 

47. The method as claimed in claim 45 wherein a time 
sequence of images frames is taken with a 80 ps gate and 
frames are taken in steps of 25 ps from -5 ns to +15 ns, and 
used for image reconstruction employing said mathematical 
inversion algorithm. 

48. The method as claimed in claim 47 wherein a suitable 
gate width between 5 ps to 1 ns and a sequence of time 
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frames from -5 to +15 ns are taken, and used for image 
reconstruction employing said mathematical inversion algo- 
rithm. 

49. The method as claimed in claim 45 wherein said 
5 mathematical inversion algorithm is modified in accordance 

with the geometry of the cell sample holder. 

50. The method as claimed in claim 45 wherein the time 
sequences of image frames are taken at angles from 0 to 90 
degrees with respect to the optical axis and used for image 
reconstruction employing the mathematical inversion algo- 

10 rithm. 

51. The method as claimed in claim 47 further comprising 
the step of detecting tumors in vivo inside human body using 
fluorescence spectroscopy. 



